Interval Truth Model

Visualizations of Main Simulation Study

Authors
Affiliation

Matthias Kloft

University of Marburg

Björn S. Siepe

University of Marburg

Daniel W. Heck

University of Marburg

Published

July 21, 2025

1 Background

This file contains the visualizations of the simulation study for the Interval Truth Model. The results of the preliminary simulation study on the two different link functions can be found in the file src/01_simulation_study_link_functions_visualizations.Qmd.

All errorbars in this document represent \(\pm 1\) Monte Carlo Standard Errors.

2 Prep

We first load all relevant packages:

Show the code
packages <- c(
  "tidyverse",
  "SimDesign",
  "here",
  "psych",
  "ggh4x",
  "ggokabeito",
  "ggExtra",
  "showtext",
  "ggdist",
  "pander",
  "sysfonts",
  "DT"
)

if (!require("pacman")) install.packages("pacman")
pacman::p_load(packages, update = F, character.only = T)

# if(!require("cmdstanr")){
#   install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
#   library(cmdstanr)
# }
  

# default chunk options
knitr::opts_chunk$set(
  fig.height = 7,
  fig.width = 10,
  include = TRUE,
  message = FALSE,
  warning = FALSE
)
source(here("src", "00_functions.R"))

# specify path to full simulation results
path_full_results <- here::here(
  "sim_results",
  "sim_results_itm_2025-04-28_03-30-27",
  "full_sim_results_itm"
)

# add google font
sysfonts::font_add_google("News Cycle", "news")
# use showtext
showtext::showtext_auto()

As the consensus locations and width have slightly different standard deviations, we can use these to standardize performance measures for ease of interpretability. Obtain the SDs from the data-generation function:

Show the code
link <- "ilr"
mean_benchmark <- simplex_to_bvn(c(.4, .2, .4), type = link)
sd_benchmark_loc <- simplex_to_bvn(c(.98, .01, .01), type = link)
sd_benchmark_wid <- simplex_to_bvn(c(.495, .01, .495), type = link)
    
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2]
    
# SD for Tr_loc
sigma_Tr_loc <- sd_benchmark_loc[1] / 4 
# SD for Tr_wid
sigma_Tr_wid <- abs(sd_benchmark_wid[2] - mean_benchmark[2]) / 4
# SD for the mean of Tr_loc and Tr_wid
sigma_Tr_interval <- sqrt(sigma_Tr_loc^2 + sigma_Tr_wid^2) / 2


# SD for a_loc on the log-scale
sigma_a_loc <- .3
# SD for b_loc
sigma_b_loc <- sigma_Tr_loc / 3
# SD for b_wid
sigma_b_wid <- sigma_Tr_wid / 3

# SD for lambda_loc and E_loc on the log-scale
sigma_lambda_E_loc <- .3
# SD for lambda_wid and E_wid on the log-scale 
sigma_lambda_E_wid <- .3

3 Deviations from Preregistration

Deviation Explanation
We used rstan instead of cmdstanr. We used rstan due to its ease of installation and use for applied users, and its useful helper functions. This should not influence the quality of our results.
We did not use multivariate priors for proficiency and discernibility parameters. TODO give reasoning.
We set \(\omega_j = 0\) in data generation. We simply forgot to mention this in the preregistration. TODO @Matze, is this correct?
Changed formula for absolute bias and MSE, with denominator changing from number of participants to the number of items.
We used 500 instead of 1,000 simulation repetitions We underestimated the computational workload, but checked the uncertainty of point estimates and realized that it was so small that 500 repetitions sufficed, which is in line with our preregistered criterion of a certain maximum MCSE.
We changed the R version. During the revision, we used the newest R version. This should not change the results.
We added additional benchmarks (i.e., simple medians). We were asked to do so during the revision, and added it because it provides an even stronger test of our model.
We added the Pearson correlation between true and estimated parameter values as an additional performance measure. We were asked to include the performance/recovery of the secondary parameter besides the consensus intervals. The correlation facilitates the interpretation of results regarding recovery.

Additionally, as we already explained in the preregistration, we did not preregister the full model setup including the priors, which left us with some freedom to choose the priors. We make these choices transparent in the manuscript.


4 Results

The full results of the simulation study are available in two data frames:

  • sim_res_itm_server.rds: Results of the simulation study, missing some MSE results due to a small coding error, but containing all information about seed, RAM usage, etc.
  • sim_res_itm_0808.rds: Resummarized results of the simulation study, now also containing all MSE results as well as more information on divergent transitions.
Show the code
sim_res_itm <- readRDS(
  here::here(
    "sim_results",
    "sim_results_itm_2025-04-28_03-30-27",
    "resummarized_results_itm.rds"
  )
)

Prepare data and convert to long format for plotting:

Show the code
sim_res_itm_long <- sim_res_itm |>
  select(!c(contains("conv"), contains("divtrans"), contains("_sd"))) |>
  # the following columns aren't needed any more for the "resummarized" results
  # "REPLICATIONS", "SIM_TIME", "RAM_USED", "SEED", "COMPLETED", "WARNINGS")) |>
  # delete "_fn_" from every column name
  rename_all( ~ str_remove(., "_fn")) |>
  pivot_longer(
    cols = !c(n_respondents, n_items),
    names_to = "measure",
    values_to = "value"
  ) |>
  # only look at non-log transformed values
  filter(!grepl("log", measure)) |>
  # rename for easier separation
  mutate(measure = str_replace(measure, "abs_bias", "absbias")) |>
  # add string so that each measure has same number of underscores
  mutate(measure = str_replace(measure, "omega", "omega_cor")) |>
  # remove only the first underscore
  mutate(measure = sub("_", "", measure, fixed = TRUE)) |>
  separate_wider_delim(measure,
                       names = c("measure", "summary", "pm", "param"),
                       delim = "_") |>
  group_by(n_respondents, n_items, measure, summary, pm) |>
  pivot_wider(names_from = "param", values_from = "value") |>
  ungroup() |>
  mutate(n_respondents = factor(n_respondents))

4.1 Convergence and Rhat

We now present the average convergence statistics. All values are averages over all replications. A performance measure with a \(>\) or a \(<\) indicates how many out of 1000 replications had a certain property. For example, \(\hat{R} < 1.1\) means that on average, a certain amount of replications had an \(\hat{R}\) value below 1.1.

Show the code
# Renaming for the table
name_mapping <- c(
  "rhat_mean" = "$\\hat{R}_{\\text{mean}}$",
  "rhat_prop1c1" = "$\\hat{R} < 1.1$",
  "rhat_prop1c05" = "$\\hat{R} < 1.05$",
  "rhat_prop1c01" = "$\\hat{R} < 1.01$",
  "ess_bulk_mean" = "$\\text{ESS}_{\\text{bulk, mean}}$",
  "ess_bulk_prop100" = "$\\text{ESS}_{\\text{bulk > 100}}$",
  "ess_bulk_prop200" = "$\\text{ESS}_{\\text{bulk > 200}}$",
  "ess_bulk_prop300" = "$\\text{ESS}_{\\text{bulk > 300}}$",
  "ess_bulk_prop400" = "$\\text{ESS}_{\\text{bulk > 400}}$",
  "ess_tail_mean" = "$\\text{ESS}_{\\text{tail, mean}}$",
  "ess_tail_prop100" = "$\\text{ESS}_{\\text{tail > 100}}$",
  "ess_tail_prop200" = "$\\text{ESS}_{\\text{tail > 200}}$",
  "ess_tail_prop300" = "$\\text{ESS}_{\\text{tail > 300}}$",
  "ess_tail_prop400" = "$\\text{ESS}_{\\text{tail > 400}}$"
)

sim_res_itm |>
  select(contains("conv")) |>
  summarize(across(everything(), mean)) |> 
  # remove "mean_conv." from every column name
  rename_all(~str_remove(., "mean_conv.")) |>
  rename_all(~str_remove(., "mcse_conv.")) |> 
  pivot_longer(cols = everything()) |> 
  # separate mean and mcse based on last underscore
  separate(name, into = c("name", "suffix"), sep = "_(?=[^_]+$)", remove = FALSE) |> 
  pivot_wider(names_from = suffix, values_from = value) |> 
  mutate(mean = round(mean, 4), 
         mcse = round(mcse, 4)) |> 
  mutate(name = name_mapping[name]) |> 
  knitr::kable()
name mean mcse
\(\hat{R}_{\text{mean}}\) 1.0011 0.0001
\(\hat{R} < 1.1\) 0.9998 0.0003
\(\hat{R} < 1.05\) 0.9997 0.0003
\(\hat{R} < 1.01\) 0.9989 0.0012
\(\text{ESS}_{\text{bulk, mean}}\) 3687.8952 23.7872
\(\text{ESS}_{\text{bulk > 100}}\) 0.9998 0.0002
\(\text{ESS}_{\text{bulk > 200}}\) 0.9998 0.0003
\(\text{ESS}_{\text{bulk > 300}}\) 0.9997 0.0004
\(\text{ESS}_{\text{bulk > 400}}\) 0.9996 0.0005
\(\text{ESS}_{\text{tail, mean}}\) 2837.8860 4.9501
\(\text{ESS}_{\text{tail > 100}}\) 1.0000 0.0001
\(\text{ESS}_{\text{tail > 200}}\) 0.9999 0.0002
\(\text{ESS}_{\text{tail > 300}}\) 0.9999 0.0002
\(\text{ESS}_{\text{tail > 400}}\) 0.9999 0.0003

The divergent transitions are shown in the table below. For purposes of readability, we omit the MCSEs here.

Show the code
sim_res_itm |>
  select(n_items, n_respondents, contains("divtrans")) |>
  select(-contains("mcse")) |>
  rename(
    "Mean DivTrans (all)" = "mean_divtrans.divergent_transitions_mean",
    "Proportion with DivTrans" = "nonz_divtrans.divergent_transitions_nonzero",
    "Mean DivTrans (models with any DivTrans)" = "nonzmean_divtrans.divergent_transitionsnonzero_mean"
  ) |>
  mutate(
    n_items = factor(n_items, ordered = TRUE),
    n_respondents = factor(n_respondents, ordered = TRUE)
  ) |>
  arrange(n_items, n_respondents) |>
  rename(Items = n_items, Respondents = n_respondents) |>
  # round numeric columns
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  knitr::kable(align = c("c", "c", "r", "r", "r"))
Items Respondents Mean DivTrans (all) Proportion with DivTrans Mean DivTrans (models with any DivTrans)
5 10 2.56 0.012 213.333
5 50 0.00 0.000 NaN
5 100 0.00 0.000 NaN
5 200 0.00 0.000 NaN
10 10 0.08 0.002 40.000
10 50 0.00 0.000 NaN
10 100 0.00 0.000 NaN
10 200 0.00 0.000 NaN
20 10 0.00 0.000 NaN
20 50 0.00 0.000 NaN
20 100 0.00 0.000 NaN
20 200 0.00 0.000 NaN
40 10 1.76 0.002 880.000
40 50 0.08 0.002 40.000
40 100 0.00 0.000 NaN
40 200 0.00 0.000 NaN

Check if our prespecified MCSE criteria (\(<.05\)) was fulfilled in all conditions:

Show the code
sim_res_itm |> 
  select(all_of(contains("mcse"))) |> 
  select(all_of(contains("bias"))) |> 
  summarize(across(everything(), max)) |> 
  pivot_longer(cols = everything()) |> 
  arrange(desc(value)) |>
  knitr::kable()
name value
E_loc_mean_fn_abs_bias_mcse 0.0217100
E_wid_mean_fn_abs_bias_mcse 0.0084091
simplemean_interval_mean_fn_abs_bias_mcse 0.0064811
E_loc_median_fn_abs_bias_mcse 0.0061761
E_wid_median_fn_abs_bias_mcse 0.0060282
simplemed_interval_mean_fn_abs_bias_mcse 0.0058701
simplemean_loc_mean_fn_abs_bias_mcse 0.0050398
Tr_interval_median_fn_abs_bias_mcse 0.0048947
Tr_interval_mean_fn_abs_bias_mcse 0.0048764
E_loc_log_mean_fn_abs_bias_mcse 0.0046608
simplemed_loc_mean_fn_abs_bias_mcse 0.0044987
lambda_wid_mean_fn_abs_bias_mcse 0.0041285
lambda_loc_mean_fn_abs_bias_mcse 0.0040510
E_loc_log_median_fn_abs_bias_mcse 0.0039721
lambda_loc_median_fn_abs_bias_mcse 0.0039011
Tr_loc_mean_fn_abs_bias_mcse 0.0038658
E_wid_log_mean_fn_abs_bias_mcse 0.0037087
simplemed_wid_mean_fn_abs_bias_mcse 0.0036033
simplemean_wid_mean_fn_abs_bias_mcse 0.0035923
Tr_loc_median_fn_abs_bias_mcse 0.0035626
lambda_wid_median_fn_abs_bias_mcse 0.0035400
lambda_wid_log_median_fn_abs_bias_mcse 0.0033802
lambda_wid_log_mean_fn_abs_bias_mcse 0.0033604
lambda_loc_log_mean_fn_abs_bias_mcse 0.0033301
lambda_loc_log_median_fn_abs_bias_mcse 0.0033097
omega_median_fn_abs_bias_mcse 0.0032760
omega_mean_fn_abs_bias_mcse 0.0031846
E_wid_log_median_fn_abs_bias_mcse 0.0030469
Tr_wid_mean_fn_abs_bias_mcse 0.0030227
a_loc_log_median_fn_abs_bias_mcse 0.0030099
a_loc_mean_fn_abs_bias_mcse 0.0029833
Tr_wid_median_fn_abs_bias_mcse 0.0029536
a_loc_log_mean_fn_abs_bias_mcse 0.0029255
a_loc_median_fn_abs_bias_mcse 0.0028233
b_loc_mean_fn_abs_bias_mcse 0.0023194
b_loc_median_fn_abs_bias_mcse 0.0023160
b_wid_mean_fn_abs_bias_mcse 0.0020543
b_wid_median_fn_abs_bias_mcse 0.0019313

4.2 Visualization: Consensus Location + Width

4.2.1 Absolute Bias of Consensus Location + Width

4.2.1.1 Combined AbsBias

Here, we show the absolute bias of the consensus location and width for the different sample sizes and number of items:

Show the code
plot_sim_absbias <- sim_res_itm_long |>
  filter(pm == "absbias") |>
  filter(summary == "mean") |>
  filter(measure %in% c("Trinterval", "simplemeaninterval", "simplemedinterval")) |>
  mutate(
    measure = case_when(
      measure == "Trinterval" ~ "ITM",
      measure == "simplemeaninterval" ~ "Simple Means",
      measure == "simplemedinterval" ~ "Simple Medians"
    )
  ) |>
  mutate(n_items = paste0(n_items, " Items")) |>
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
  # # Standardize with true variability
  # mutate(mean = mean / sigma_Tr_interval,
  #        mcse = mcse / sigma_Tr_interval) |>
  ggplot(aes(
    x = n_respondents,
    y = mean,
    color = measure,
    group = measure
  )) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
  geom_point(position = position_dodge(0.7), size = 2.5) +
  geom_errorbar(
    aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
    width = 0.95,
    position = position_dodge(0.7),
    show.legend = FALSE
  ) +
  ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
  scale_y_continuous(limits = c(0, .49), expand = c(0, 0)) +
  ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
  labs(x = "Number of Respondents", y = "Absolute Bias", color = "") +
  theme_icm() +
  theme(legend.position = "top", text = element_text(size = 22))
# For manuscript plot:
# theme(legend.position = "top",
#       text = element_text(size = 14),
#       axis.text.x = element_text(size = 13, margin = margin(b = 10)),
#       axis.text.y = element_text(size = 13, margin = margin(l = 10)),
#       legend.text = element_text(size = 15.5))
#
ggsave(
  here("plots", "sim_main", "sim_absbias.pdf"),
  plot_sim_absbias,
  width = 9,
  height = 4.5
)

plot_sim_absbias

4.2.1.2 Separate AbsBias

Here, we focus on the ITM and separate the bias by location and width, while standardizing for true variability:

Show the code
plot_sim_bias_widloc <- sim_res_itm_long |>
  filter(pm == "absbias") |>
  filter(summary == "mean") |> 
  filter(measure %in% c("Trloc", "Trwid")) |> 
  mutate(measure = case_when(
    measure == "Trloc" ~ "Location",
    measure == "Trwid" ~ "Width"
  )) |>
  mutate(n_items = paste0(n_items, " Items")) |> 
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", 
                                              "10 Items", 
                                              "20 Items",
                                              "40 Items"))) |>
  # Standardize with true variability
  mutate(mean = if_else(measure == "Location",
                        mean / sigma_Tr_loc,
                        mean / sigma_Tr_wid
                        ),
         mcse = if_else(measure == "Location",
                        mcse / sigma_Tr_loc,
                        mcse / sigma_Tr_wid
                        )) |>
  ggplot(aes(x = n_respondents, 
             y = mean, 
             color = measure,
             group = measure)) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1))+
  geom_point(position = position_dodge(0.7), 
             size = 2.5) +
  geom_errorbar(aes(ymin = mean - 1*mcse,
                            ymax = mean + 1*mcse),
                        width = .8,
                 position = position_dodge(0.7),
                 show.legend = FALSE)+
  ggh4x::facet_wrap2(n_items ~ .,
                     axes = "all",
                     nrow = 1) +
  scale_y_continuous(limits = c(0, 0.31), expand = c(0,0)) +
  ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
  labs(x = "Number of Respondents",
       y = "Standardized Absolute Bias",
       color = "") +
  theme_icm()+
  theme(legend.position = "top",
        text = element_text(size = 28))

ggsave(here("plots","sim_main", "sim_absbias_widloc.pdf"),
       plot_sim_bias_widloc, width = 7, height = 3.75)

plot_sim_bias_widloc

Alternatively, combine plots for location and width and standardize with true variability:

Show the code
plot_sim_absbias_widloc_comb <- sim_res_itm_long |>
  filter(pm == "absbias") |>
  filter(summary == "mean") |>
  filter(
    measure %in% c(
      "Trloc",
      "Trwid",
      "simplemeanloc",
      "simplemeanwid",
      "simplemedloc",
      "simplemedwid"
    )
  ) |>
  # split based on model
  mutate(
    model = case_when(
      grepl("^Tr", measure) ~ "ITM",
      grepl("^simplemean", measure) ~ "Simple Means",
      grepl("^simplemed", measure) ~ "Simple Medians"
    ),
    measure = sub("^Tr|simplemean|simplemed", "", measure)
  ) |>
  mutate(measure = case_when(measure == "loc" ~ "Location", measure == "wid" ~ "Width")) |>
  # Standardize with true variability
  mutate(
    mean = if_else(measure == "Location", mean / sigma_Tr_loc, mean / sigma_Tr_wid),
    mcse = if_else(measure == "Location", mcse / sigma_Tr_loc, mcse / sigma_Tr_wid)
  ) |>
  mutate(n_items = paste0(n_items, " Items")) |>
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
  ggplot(aes(
    x = n_respondents,
    y = mean,
    color = model,
    group = model
  )) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
  geom_point(position = position_dodge(0.7), size = 2.5) +
  geom_errorbar(
    aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
    width = .8,
    position = position_dodge(0.7),
    show.legend = FALSE
  ) +
  ggh4x::facet_grid2(measure ~ n_items, axes = "all") +
  scale_y_continuous(limits = c(0, 0.39), expand = c(0, 0)) +
  ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
  labs(x = "Number of Respondents", y = "Standardized Absolute Bias", color = "") +
  theme_icm() +
  # theme(legend.position = "top", text = element_text(size = 14))+
  # For manuscript plot:
  theme(
    legend.position = "top",
    text = element_text(size = 15.5),
    axis.text.x = element_text(size = 14., margin = margin(b = 10)),
    axis.text.y = element_text(size = 14.5, margin = margin(l = 10)),
    legend.text = element_text(size = 17.5)
  )

plot_sim_absbias_widloc_comb

Show the code
ggsave(here("plots","sim_main", "sim_absbias_widloc_comb.pdf"),
       plot_sim_absbias_widloc_comb, width = 11, height = 7)

4.2.2 MSE for Location and Width

4.2.2.1 Combined MSE

Show both ITM and simple means/medians in a comparison:

Show the code
plot_sim_mse <- sim_res_itm_long |>
  filter(pm == "mse") |>
  filter(summary == "mean") |> 
  filter(measure %in% c("Trinterval", "simplemeaninterval", "simplemedinterval")) |> 
  mutate(measure = case_when(
    measure == "Trinterval" ~ "ITM",
    measure == "simplemeaninterval" ~ "Simple Means",
    measure == "simplemedinterval" ~ "Simple Medians"
  )) |>
  mutate(n_items = paste0(n_items, " Items")) |> 
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", 
                                              "10 Items", 
                                              "20 Items",
                                              "40 Items"))) |>
  ggplot(aes(x = n_respondents, 
             y = mean, 
             color = measure,
             group = measure)) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1))+
  geom_point(position = position_dodge(0.7), 
             size = 2.5) +
  geom_errorbar(aes(ymin = mean - 1*mcse,
                            ymax = mean + 1*mcse),
                        width = .8,
                 position = position_dodge(0.7),
                 show.legend = FALSE)+
  ggh4x::facet_wrap2(n_items ~ .,
                     axes = "all",
                     nrow = 1) +
  scale_y_continuous(limits = c(0, .49), expand = c(0,0)) +
  ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3))+
  labs(x = "Number of Respondents",
       y = "MSE",
       color = "") +
  theme_icm()+
  theme(legend.position = "top",
        text = element_text(size = 28))

ggsave(here("plots","sim_main", "sim_mse.pdf"),
       plot_sim_mse, width = 7, height = 3.75)

plot_sim_mse

4.2.2.2 Separate MSE

Separated by location and width and standardize values with the true variance:

Show the code
plot_sim_mse_widloc <- sim_res_itm_long |>
  filter(pm == "mse") |>
  filter(summary == "mean") |>
  filter(measure %in% c("Trloc", "Trwid")) |>
  mutate(measure = case_when(measure == "Trloc" ~ "Location", measure == "Trwid" ~ "Width")) |>
  mutate(n_items = paste0(n_items, " Items")) |>
  # standardize with true variability
  mutate(
    mean = if_else(
      measure == "Location",
      mean / sigma_Tr_loc^2,
      mean / sigma_Tr_wid^2
    ),
    mcse = if_else(
      measure == "Location",
      mcse / sigma_Tr_loc^2,
      mcse / sigma_Tr_wid^2
    )
  ) |>
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
  ggplot(aes(
    x = n_respondents,
    y = mean,
    color = measure,
    group = measure
  )) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
  geom_point(position = position_dodge(0.7), size = 2.5) +
  geom_errorbar(
    aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
    width = .8,
    position = position_dodge(0.7),
    show.legend = FALSE
  ) +
  ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
  scale_y_continuous(limits = c(0, 0.151), expand = c(0, 0)) +
  ggokabeito::scale_color_okabe_ito(order = c(2, 7)) +
  labs(x = "Number of Respondents", y = "MSE", color = "") +
  theme_icm() +
  theme(legend.position = "top", text = element_text(size = 28))

ggsave(here("plots","sim_main", "sim_mse_widloc.pdf"),
       plot_sim_mse_widloc, width = 7, height = 3.75)

plot_sim_mse_widloc

Alternatively, combine plots for location and width and standardize values with the true variance:

Show the code
plot_sim_mse_widloc_comb <- sim_res_itm_long |>
  filter(pm == "mse") |>
  filter(summary == "mean") |>
  filter(measure %in% c("Trloc", "Trwid", "simplemeanloc", "simplemeanwid", "simplemedloc", "simplemedwid")) |>
  # split based on model
  mutate(model = case_when(
    grepl("^Tr", measure) ~ "ITM",
    grepl("^simplemean", measure) ~ "Simple Means",
    grepl("^simplemed", measure) ~ "Simple Medians"
  ),
  measure = sub("^Tr|simplemean|simplemed", "", measure)) |> 
  mutate(measure = case_when(
    measure == "loc" ~ "Location", 
    measure == "wid" ~ "Width")) |>
  mutate(n_items = paste0(n_items, " Items")) |>
  # standardize with true variability
  mutate(mean = if_else(measure == "Location",
                        mean / sigma_Tr_loc^2,
                        mean / sigma_Tr_wid^2
                        ),
         mcse = if_else(measure == "Location",
                        mcse / sigma_Tr_loc^2,
                        mcse / sigma_Tr_wid^2
                        )) |>
  # order n_items correctly
  mutate(n_items = factor(n_items, 
                          levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
  ggplot(aes(
    x = n_respondents,
    y = mean,
    color = model,
    group = model
  )) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
  geom_point(position = position_dodge(0.7), size = 2.5) +
  geom_errorbar(
    aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
    width = .8,
    position = position_dodge(0.7),
    show.legend = FALSE
  ) +
  ggh4x::facet_grid2(measure ~ n_items, 
                     axes = "all") +
  scale_y_continuous(limits = c(0, 0.21), expand = c(0, 0)) +
  ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
  labs(x = "Number of Respondents", y = "Standardized MSE", color = "") +
  theme_icm() +
  # theme(legend.position = "top", text = element_text(size = 14))+
  # For manuscript plot:
  theme(legend.position = "top",
        text = element_text(size = 14),
        axis.text.x = element_text(size = 13.5, margin = margin(b = 10)),
        axis.text.y = element_text(size = 13.5, margin = margin(l = 10)),
        legend.text = element_text(size = 15.5))

plot_sim_mse_widloc_comb

Show the code
ggsave(here("plots","sim_main", "sim_mse_widloc_comb.pdf"),
       plot_sim_mse_widloc_comb, width = 11, height = 7)

4.2.3 Scatterplot of Location and Width Bias

Apply a function that retrieves location and width estimates to the results:

Show the code
if(exists(here(
  "sim_results",
  "sim_results_itm_2025-04-28_03-30-27",
  "locwid_bias.rds"
))) {
  # read the data
  locwid_bias <- readRDS(here(
    "sim_results",
    "sim_results_itm_2025-04-28_03-30-27",
    "locwid_bias.rds"
  ))
} else{
  locwid_bias <- prep_locwid(path_full_results)
  # save the data
  saveRDS(
    locwid_bias,
    here(
      "sim_results",
      "sim_results_itm_2025-04-28_03-30-27",
      "locwid_bias.rds"
    )
  )
}

Then create a scatter plot to investigate compensatory behavior, in other words, if the bias of the location is higher when the bias of the width is lower and vice versa. Overall, there does not seem to be strong evidence for compensatory behavior.

Show the code
#  standardize the values
locwid_bias <- locwid_bias |>
  dplyr::mutate(loc_bias_std = loc_bias / sigma_Tr_loc,
                wid_bias_std = wid_bias / sigma_Tr_wid)


# if fonts aren't working properly with ggExtra
# extrafont::font_import()

# # add google font
sysfonts::font_add_google("News Cycle", "news")
# use showtext
showtext::showtext_auto()
# windows()
# extrafont::font_import(pattern = "NewsCycle-Regular.ttf")
# extrafont::loadfonts()



scatter_tmp <- locwid_bias |>
  ggplot(aes(x = loc_bias_std, y = wid_bias_std)) +
  geom_point(alpha = 0.5) +
  geom_abline(intercept = 0,
              slope = 1,
              linetype = "dashed") +
  labs(x = "Location Bias (standardized)", y = "Width Bias (standardized)") +
  geom_smooth(method = "loess", se = TRUE) +
  scale_x_continuous(limits = c(0, .85), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, .85), expand = c(0, 0)) +
  # manually set theme here due to font conflicts
  ggplot2::theme_minimal(base_family = "news") +
  ggplot2::theme(
    # remove minor grid
    panel.grid.minor = ggplot2::element_blank(),
    # Title and Axis Texts
    plot.title = ggplot2::element_text(
      face = "plain",
      size = ggplot2::rel(1.2),
      hjust = 0.5
    ),
    plot.subtitle = ggplot2::element_text(size = ggplot2::rel(1.1), hjust = 0.5),
    axis.text.x = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.05)),
    axis.text.y = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.05)),
    axis.title.x = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.3)),
    axis.title.y = ggplot2::element_text(face = "plain", size = ggplot2::rel(1.3)),
    axis.line = element_line(colour = "#6d6d6e"),
    
    # Faceting
    strip.text = ggplot2::element_text(
      face = "plain",
      size = ggplot2::rel(1.1),
      hjust = 0.5
    ),
    strip.text.x.top = ggplot2::element_text(
      face = "plain",
      size = ggplot2::rel(1.2),
      hjust = 0.5
    ),
    # strip.text.y = element_blank(),
    strip.background = ggplot2::element_rect(fill = NA, color = NA),
    # Grid
    panel.grid = ggplot2::element_line(colour = "#F3F4F5"),
    # Legend
    legend.title = ggplot2::element_text(face = "plain"),
    legend.position = "top",
    legend.justification = 1,
    # Panel/Facets
    panel.spacing.x = ggplot2::unit(1.6, "lines"),
    panel.spacing.y = ggplot2::unit(1.6, "lines"),
    # Remove vertical grid lines
    panel.grid.major.x = ggplot2::element_blank()
  ) +
  theme(text = element_text(size = 28))

# show marginal distributions as histograms
# doesn't work with custom font
# scatter_locwid <- ggExtra::ggMarginal(scatter_tmp, type = "histogram", bins = 50)+
#   theme_minimal(base_family = "news")

ggsave(here("plots","sim_main", "scatter_locwid.pdf"),
       scatter_tmp, width = 7, height = 5)

scatter_tmp

4.2.4 Visualization of Raw Results

Again, apply a function that retrieves location and width estimates to the results:

Show the code
if(exists(here(
  "sim_results",
  "sim_results_itm_2025-04-28_03-30-27",
  "raw_absbias.rds"
))) {
  # read the data
  path_full_results <- here("sim_results",
                            "sim_results_itm_2025-04-28_03-30-27",
                            "raw_absbias.rds")
} else{
  # function to obtain the raw results
  raw_absbias <-
    prep_locwid(path_full_results,
                simplemeans = TRUE,
                simplemedians = TRUE)
  # save the data
  saveRDS(
    raw_absbias,
    here(
      "sim_results",
      "sim_results_itm_2025-04-28_03-30-27",
      "raw_absbias.rds"
    )
  )
}

Or just read in the results:

Show the code
raw_absbias <- readRDS(here("sim_results",
                            "sim_results_itm_2025-04-28_03-30-27",
                            "raw_absbias.rds"))

Instead of only showing aggregate performance measures, we can also visualize the raw repetition-wise performance measures to investigate the variability and potential skewness of the performance.

Show the code
# attach them to the condition-wise results
raw_ests_boxplot <- sim_res_itm |>
  mutate(iteration = row_number()) |>
  dplyr::select(n_respondents, n_items, iteration) |>
  right_join(raw_absbias, by = "iteration") |>
  pivot_longer(cols = c(loc_bias, wid_bias, smloc_bias, smwid_bias, smmloc_bias, smmwid_bias)) |>
  mutate(n_respondents = factor(n_respondents)) |>
  mutate(name = gsub("_bias", "", name)) |>
  mutate(
    name = case_match(
      name,
      "loc" ~ "ITM Location",
      "wid" ~ "ITM Width",
      "smloc" ~ "Simplemeans Location",
      "smwid" ~ "Simplemeans Width",
      "smmloc" ~ "Simplemedians Location",
      "smmwid" ~ "Simplemedians Width",
      .default = name
    )
  ) |>
  separate_wider_delim(name, delim = " ", names = c("model", "measure")) |>
  mutate(model = gsub("Simplemeans", "Simple Means", model)) |>
  mutate(model = gsub("Simplemedians", "Simple Medians", model)) |>
  mutate(n_items = paste0(n_items, " Items")) |>
  # order n_items correctly
  mutate(n_items = factor(n_items, levels = c("5 Items", "10 Items", "20 Items", "40 Items"))) |>
  mutate(n_respondents = as_factor(n_respondents)) |>
  ggplot(aes(
    x = n_respondents,
    y = value,
    color = model,
    group = model,
    fill = model
  )) +
  # add vertical line between different sample sizes
  geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
  # ggdist::stat_halfeye(
  #   position = "dodge",
  #   scale = 5
  # ) +
  geom_boxplot(
    alpha = .2,
    lwd = 0.6,
    fatten = 0.9,
    inherit.aes = FALSE,
    aes(
      x = n_respondents,
      y = value,
      fill = model,
      color = model
    )
  ) +
  # geom_point(
  #   size = 1,
  #   alpha = 0.15,
  #   position = ggplot2::position_dodge(width = .6)
  # )+
  # geom_jitter(size=0.4, alpha=0.2) +
  ggh4x::facet_grid2(measure ~ n_items, axes = "all", scales = "free_y") +
  # scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
  ggokabeito::scale_color_okabe_ito(order = c(5, 1, 3)) +
  ggokabeito::scale_fill_okabe_ito(order = c(5, 1, 3)) +
  labs(x = "Number of Respondents",
       y = "Absolute Bias",
       color = "",
       fill = "") +
  theme_icm() +
  # theme(legend.position = "top", text = element_text(size = 18))+
  # For manuscript plot:
  theme(
    legend.position = "top",
    text = element_text(size = 14),
    axis.text.x = element_text(size = 13.5, margin = margin(b = 10)),
    axis.text.y = element_text(size = 13.5, margin = margin(l = 10)),
    legend.text = element_text(size = 15.5)
  )

raw_ests_boxplot

Show the code
ggsave(here("plots","sim_main", "raw_ests_boxplot.pdf"),
       raw_ests_boxplot, width = 15, height = 7)

4.2.5 Summary Table

To show the numerical results of the main outcome measures, we create a summary table below. It shows unstandardized estimates.

Show the code
sim_res_itm |>
  select(
    n_items,
    n_respondents,
    Tr_interval_mean_fn_abs_bias_mean,
    Tr_interval_mean_fn_abs_bias_mcse,
    Tr_interval_mean_fn_mse_mean,
    Tr_interval_mean_fn_mse_mcse,
    simplemean_interval_mean_fn_abs_bias_mcse,
    simplemean_interval_mean_fn_abs_bias_mean,
    simplemean_interval_mean_fn_mse_mean,
    simplemean_interval_mean_fn_mse_mcse,
    simplemed_interval_mean_fn_abs_bias_mcse,
    simplemed_interval_mean_fn_abs_bias_mean,
    simplemed_interval_mean_fn_mse_mean,
    simplemed_interval_mean_fn_mse_mcse
  ) |>
  pivot_longer(cols = !c(n_items, n_respondents)) |>
  # split based on last underscore
  separate(
    name,
    into = c("name", "suffix"),
    sep = "_(?=[^_]+$)",
    remove = FALSE
  ) |>
  pivot_wider(names_from = suffix, values_from = value) |>
  # again split based on last underscore
  separate(name,
           into = c("name", "pm"),
           sep = "_(?=[^_]+$)",
           remove = FALSE) |>
  # remove "mean_fn" from name
  mutate(name = str_remove(name, "_mean_fn_abs")) |>
  mutate(name = str_remove(name, "_mean_fn")) |>
  dplyr::rename(
    "Items" = "n_items",
    "Respondents" = "n_respondents",
    "Model" = "name",
    "Measure" = "pm",
    "Mean" = "mean",
    "MCSE" = "mcse"
  ) |>
  mutate(
    Model = case_when(
      Model == "Tr_interval" ~ "ITM",
      Model == "simplemean_interval" ~ "Simple Means",
      Model == "simplemed_interval" ~ "Simple Medians"
    )
  ) |>
  mutate(Measure = case_when(Measure == "bias" ~ "Bias", Measure == "mse" ~ "MSE")) |>
  arrange(Items, Respondents) |>
  mutate(across(c(Mean, MCSE), ~ round(., 4))) |>
  knitr::kable(align = c("c", "c", "l", "l", "c", "c"))
Items Respondents Model Measure Mean MCSE
5 10 ITM Bias 0.4117 0.0049
5 10 ITM MSE 0.1445 0.0035
5 10 Simple Means Bias 0.4708 0.0065
5 10 Simple Means MSE 0.1995 0.0060
5 10 Simple Medians Bias 0.4838 0.0059
5 10 Simple Medians MSE 0.4838 0.0061
5 50 ITM Bias 0.1870 0.0024
5 50 ITM MSE 0.0302 0.0008
5 50 Simple Means Bias 0.2114 0.0027
5 50 Simple Means MSE 0.0390 0.0010
5 50 Simple Medians Bias 0.2294 0.0029
5 50 Simple Medians MSE 0.2294 0.0030
5 100 ITM Bias 0.1323 0.0016
5 100 ITM MSE 0.0149 0.0004
5 100 Simple Means Bias 0.1541 0.0020
5 100 Simple Means MSE 0.0205 0.0005
5 100 Simple Medians Bias 0.1615 0.0019
5 100 Simple Medians MSE 0.1615 0.0020
5 200 ITM Bias 0.0941 0.0012
5 200 ITM MSE 0.0076 0.0002
5 200 Simple Means Bias 0.1127 0.0015
5 200 Simple Means MSE 0.0112 0.0003
5 200 Simple Medians Bias 0.1129 0.0014
5 200 Simple Medians MSE 0.1129 0.0014
10 10 ITM Bias 0.3945 0.0037
10 10 ITM MSE 0.1321 0.0026
10 10 Simple Means Bias 0.4655 0.0046
10 10 Simple Means MSE 0.1908 0.0043
10 10 Simple Medians Bias 0.4780 0.0047
10 10 Simple Medians MSE 0.4780 0.0046
10 50 ITM Bias 0.1801 0.0015
10 50 ITM MSE 0.0276 0.0005
10 50 Simple Means Bias 0.2139 0.0020
10 50 Simple Means MSE 0.0404 0.0009
10 50 Simple Medians Bias 0.2231 0.0019
10 50 Simple Medians MSE 0.2231 0.0019
10 100 ITM Bias 0.1296 0.0012
10 100 ITM MSE 0.0145 0.0003
10 100 Simple Means Bias 0.1553 0.0015
10 100 Simple Means MSE 0.0210 0.0004
10 100 Simple Medians Bias 0.1618 0.0014
10 100 Simple Medians MSE 0.1618 0.0015
10 200 ITM Bias 0.0908 0.0008
10 200 ITM MSE 0.0070 0.0001
10 200 Simple Means Bias 0.1124 0.0011
10 200 Simple Means MSE 0.0110 0.0002
10 200 Simple Medians Bias 0.1144 0.0010
10 200 Simple Medians MSE 0.1144 0.0010
20 10 ITM Bias 0.3934 0.0029
20 10 ITM MSE 0.1326 0.0019
20 10 Simple Means Bias 0.4618 0.0035
20 10 Simple Means MSE 0.1872 0.0031
20 10 Simple Medians Bias 0.4744 0.0034
20 10 Simple Medians MSE 0.4744 0.0033
20 50 ITM Bias 0.1750 0.0013
20 50 ITM MSE 0.0263 0.0004
20 50 Simple Means Bias 0.2129 0.0014
20 50 Simple Means MSE 0.0397 0.0006
20 50 Simple Medians Bias 0.2220 0.0014
20 50 Simple Medians MSE 0.2220 0.0014
20 100 ITM Bias 0.1266 0.0008
20 100 ITM MSE 0.0136 0.0002
20 100 Simple Means Bias 0.1555 0.0011
20 100 Simple Means MSE 0.0211 0.0003
20 100 Simple Medians Bias 0.1605 0.0010
20 100 Simple Medians MSE 0.1605 0.0010
20 200 ITM Bias 0.0892 0.0006
20 200 ITM MSE 0.0067 0.0001
20 200 Simple Means Bias 0.1131 0.0008
20 200 Simple Means MSE 0.0112 0.0002
20 200 Simple Medians Bias 0.1141 0.0007
20 200 Simple Medians MSE 0.1141 0.0007
40 10 ITM Bias 0.4018 0.0022
40 10 ITM MSE 0.1375 0.0016
40 10 Simple Means Bias 0.4625 0.0030
40 10 Simple Means MSE 0.1885 0.0027
40 10 Simple Medians Bias 0.4757 0.0026
40 10 Simple Medians MSE 0.4757 0.0026
40 50 ITM Bias 0.1754 0.0009
40 50 ITM MSE 0.0262 0.0003
40 50 Simple Means Bias 0.2124 0.0012
40 50 Simple Means MSE 0.0390 0.0004
40 50 Simple Medians Bias 0.2233 0.0011
40 50 Simple Medians MSE 0.2233 0.0011
40 100 ITM Bias 0.1233 0.0006
40 100 ITM MSE 0.0129 0.0001
40 100 Simple Means Bias 0.1531 0.0008
40 100 Simple Means MSE 0.0204 0.0002
40 100 Simple Medians Bias 0.1586 0.0007
40 100 Simple Medians MSE 0.1586 0.0008
40 200 ITM Bias 0.0865 0.0004
40 200 ITM MSE 0.0064 0.0001
40 200 Simple Means Bias 0.1117 0.0007
40 200 Simple Means MSE 0.0110 0.0001
40 200 Simple Medians Bias 0.1135 0.0005
40 200 Simple Medians MSE 0.1135 0.0005

4.3 Absolute Bias and MSE for Other Parameters

The performance measures of all other parameters can be accessed in the tabs below.

Create a plotting function:

Show the code
plot_bias <- function(data, pm, param, lab_y_axis = "Absolute Bias", ...) {
  # obtain y-axis scaling
  y_max <- data |>
    filter(pm == {{pm}}) |>
    filter(summary == "mean") |>
    filter(measure == {{param}}) |>
    mutate(max_val = mean + 1 * mcse) |>
    pull(max_val)
  y_max <- max(y_max) * 1.1
  
  data |>
    dplyr::filter(pm == {{pm}}) |>
    dplyr::filter(summary == "mean") |>
    dplyr::filter(measure == {{param}}) |>
    dplyr::mutate(
      measure = case_when(
        measure == "Trloc" ~ "Consensus Location",
        measure == "Trwid" ~ "Consensus Width",
        measure == "Trinterval" ~ "Interval",
        measure == "aloc" ~ "a Location",
        measure == "bloc" ~ "b Location",
        measure == "bwid" ~ "b Width",
        measure == "Eloc" ~ "E Location",
        measure == "Ewid" ~ "E Width",
        measure == "lambdaloc" ~ "Lambda Location",
        measure == "lambdawid" ~ "Lambda Width",
        measure == "omegacor" ~ "Omega Correlation",
        measure == "simplmeanloc" ~ "Simplemeans Location",
        measure == "simplmeanwid" ~ "Simplemeans Width",
        measure == "simplemeaninterval" ~ "Simplemeans Interval"
      )
    ) |>
    dplyr::mutate(n_items = paste0(n_items, " Items")) |>
    # order n_items correctly
    dplyr::mutate(n_items = factor(n_items, levels = c(
      "5 Items", "10 Items", "20 Items", "40 Items"
    ))) |>
    ggplot(aes(x = n_respondents, y = mean)) +
    # add vertical line between different sample sizes
    geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
    geom_point(size = 2.5) +
    geom_errorbar(
      aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
      width = .8,
      show.legend = FALSE
    ) +
    ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
    scale_y_continuous(expand = c(0, 0), limits = c(0, y_max)) +
    # ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
    labs(x = "Number of Respondents", y = lab_y_axis) +
    theme_icm() +
    theme(legend.position = "top", text = element_text(size = 28))
}

4.3.1 Absolute Bias

Additionally, we can plot the standardized absolute biases and MSEs on the log scale.

First, prep the data to select parameters on the log-scale:

Show the code
sim_res_log <- sim_res_itm |>
  select(
    n_respondents,
    n_items,
    contains("log"),
    contains("omega"),
    contains("Tr_loc"),
    contains("Tr_wid"),
    contains("b_loc"),
    contains("b_wid")
  ) |>
  rename_all(~ str_remove(., "_fn")) |>
  pivot_longer(
    cols = !c(n_respondents, n_items),
    names_to = "measure",
    values_to = "value"
  ) |>
  # remove the "log", as everything is on a log scale (besides omegacor)
  mutate(measure = str_remove(measure, "log_")) |>
  # rename for easier separation
  mutate(measure = str_replace(measure, "abs_bias", "absbias")) |>
  # add string so that each measure has same number of underscores
  mutate(measure = str_replace(measure, "omega", "omega_cor")) |>
  # remove only the first underscore
  mutate(measure = sub("_", "", measure, fixed = TRUE)) |>
  separate_wider_delim(measure,
                       names = c("measure", "summary", "pm", "param"),
                       delim = "_") |>
  group_by(n_respondents, n_items, measure, summary, pm) |>
  pivot_wider(names_from = "param", values_from = "value") |>
  ungroup() |>
  mutate(n_respondents = factor(n_respondents))

4.3.2 MSE

4.4 Correlation of Estimates with True Parameters

Define the plotting function for correlations

Show the code
plot_corr <- function(data, pm, param, ...) {
  # obtain y-axis scaling
  y_max <- data |>
    filter(pm == {{pm}}) |>
    filter(summary == "mean") |>
    filter(measure == {{param}}) |>
    mutate(max_val = mean + 1 * mcse) |>
    pull(max_val)
  y_max <- max(y_max) + .005
  
  y_min <- data |>
    filter(pm == {{pm}}) |>
    filter(summary == "mean") |>
    filter(measure == {{param}}) |>
    mutate(min_val = mean - 1 * mcse) |>
    pull(min_val)
  y_min <- min(y_min) - .005
  
  # plot
  data |>
    dplyr::filter(pm == {{pm}}) |>
    dplyr::filter(summary == "mean") |>
    dplyr::filter(measure == {{param}}) |>
    dplyr::mutate(
      measure = case_when(
        measure == "Trloc" ~ "Consensus Location",
        measure == "Trwid" ~ "Consensus Width",
        measure == "aloc" ~ "a Location (log-scale)",
        measure == "bloc" ~ "b Location",
        measure == "bwid" ~ "b Width",
        measure == "Eloc" ~ "E Location (log-scale)",
        measure == "Ewid" ~ "E Width (log-scale)",
        measure == "lambdaloc" ~ "Lambda Location (log-scale)",
        measure == "lambdawid" ~ "Lambda Width (log-scale)",
        measure == "omegacor" ~ "Omega Correlation"
      )
    ) |>
    dplyr::mutate(n_items = paste0(n_items, " Items")) |>
    # order n_items correctly
    dplyr::mutate(n_items = factor(n_items, levels = c(
      "5 Items", "10 Items", "20 Items", "40 Items"
    ))) |>
    ggplot(aes(x = n_respondents, y = mean)) +
    # add vertical line between different sample sizes
    geom_vline(colour = "#F3F4F5", xintercept = seq(1.5, 4, 1)) +
    geom_hline(yintercept = 1,
               linetype = "dashed",
               color = "grey50") +
    geom_point(size = 2.5) +
    geom_errorbar(
      aes(ymin = mean - 1 * mcse, ymax = mean + 1 * mcse),
      width = .8,
      show.legend = FALSE
    ) +
    ggh4x::facet_wrap2(n_items ~ ., axes = "all", nrow = 1) +
    scale_y_continuous(expand = expansion(), limits = c(y_min, y_max)) +
    # ggokabeito::scale_color_okabe_ito(order = c(2, 7))+
    labs(x = "Number of Respondents", y = "Absolute Bias") +
    theme_icm() +
    theme(legend.position = "top", text = element_text(size = 28))
}

4.5 All Measures for All Parameters in One Summary Table

Show the code
sim_res_log |>
  dplyr::filter(measure %in% rel_params,
                summary == "mean",
                pm %in% c("absbias", "mse", "corr")) |>
  # round mean and mcse
  dplyr::mutate(mean = round(mean, 3), mcse = round(mcse, 3)) |>
  # pivot_wider of "pm"
  pivot_wider(names_from = pm, values_from = c(mean, mcse, sd)) |>
  # reorder columns
  select(
    n_respondents,
    n_items,
    measure,
    mean_absbias,
    mcse_absbias,
    mean_mse,
    mcse_mse,
    mean_corr,
    mcse_corr
  ) |>
  # replace values in measure with proper names
  mutate(
    measure = case_when(
      measure == "Trloc" ~ "Consensus Location",
      measure == "Trwid" ~ "Consensus Width",
      measure == "aloc" ~ "a Location (log-scale)",
      measure == "bloc" ~ "b Location",
      measure == "bwid" ~ "b Width",
      measure == "Eloc" ~ "E Location (log-scale)",
      measure == "Ewid" ~ "E Width (log-scale)",
      measure == "lambdaloc" ~ "Lambda Location (log-scale)",
      measure == "lambdawid" ~ "Lambda Width (log-scale)",
      measure == "omegacor" ~ "Omega Correlation",
      TRUE ~ measure
    )
  ) |>
  rename(
    "Items" = n_items,
    "Respondents" = n_respondents,
    Parameter = measure,
    "Mean Absolute Bias" = mean_absbias,
    "MCSE Absolute Bias" = mcse_absbias,
    "Mean MSE" = mean_mse,
    "MCSE MSE" = mcse_mse,
    "Mean Correlation" = mean_corr,
    "MCSE Correlation" = mcse_corr
  ) |>
  mutate(
    Items = factor(Items),
    Respondents = factor(Respondents),
    Parameter = factor(Parameter)
  ) |>
  arrange(Items, Respondents) |>
  #knitr::kable(align = c("c", "c", "l", "c", "c", "c", "c"))
  DT::datatable(
    filter = "top",
    rownames = FALSE,
    options = list(
      pageLength = 50,
      lengthMenu = c(50, 100),
      dom = "t",
      columnDefs = list(list(
        className = 'dt-center',
        targets = c(
          "Items",
          "Respondents",
          "Mean Absolute Bias",
          "MCSE Absolute Bias",
          "Mean MSE",
          "MCSE MSE",
          "Mean Correlation",
          "MCSE Correlation"
        )
      )),
      autoWidth = TRUE,
      scrollX = TRUE,
      scrollY = "720px"
    )
  )

5 Session Info

Show the code
pander::pander(sessionInfo())

R version 4.5.0 (2025-04-11 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=German_Germany.utf8, LC_CTYPE=German_Germany.utf8, LC_MONETARY=German_Germany.utf8, LC_NUMERIC=C and LC_TIME=German_Germany.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: extraDistr(v.1.10.0), DT(v.0.33), pander(v.0.6.6), ggdist(v.3.3.3), showtext(v.0.9-7), showtextdb(v.3.0), sysfonts(v.0.8.9), ggExtra(v.0.10.1), ggokabeito(v.0.1.0), ggh4x(v.0.3.1), psych(v.2.5.6), here(v.1.0.1), SimDesign(v.2.19.2), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.3.0), ggplot2(v.3.5.2), tidyverse(v.2.0.0) and pacman(v.0.5.1)

loaded via a namespace (and not attached): tidyselect(v.1.2.1), farver(v.2.1.2), R.utils(v.2.13.0), fastmap(v.1.2.0), promises(v.1.3.3), digest(v.0.6.37), timechange(v.0.3.0), mime(v.0.13), lifecycle(v.1.0.4), magrittr(v.2.0.3), compiler(v.4.5.0), sass(v.0.4.10), rlang(v.1.1.6), tools(v.4.5.0), yaml(v.2.3.10), knitr(v.1.50), labeling(v.0.4.3), htmlwidgets(v.1.6.4), mnormt(v.2.1.1), curl(v.6.4.0), RColorBrewer(v.1.1-3), miniUI(v.0.1.2), withr(v.3.0.2), R.oo(v.1.27.1), grid(v.4.5.0), xtable(v.1.8-4), future(v.1.58.0), progressr(v.0.15.1), globals(v.0.18.0), scales(v.1.4.0), cli(v.3.6.5), rmarkdown(v.2.29), ragg(v.1.4.0), generics(v.0.1.4), rstudioapi(v.0.17.1), future.apply(v.1.20.0), tzdb(v.0.5.0), sessioninfo(v.1.2.3), cachem(v.1.1.0), pbapply(v.1.7-2), splines(v.4.5.0), audio(v.0.1-11), parallel(v.4.5.0), beepr(v.2.0), vctrs(v.0.6.5), Matrix(v.1.7-3), jsonlite(v.2.0.0), hms(v.1.1.3), listenv(v.0.9.1), crosstalk(v.1.2.1), systemfonts(v.1.2.3), testthat(v.3.2.3), jquerylib(v.0.1.4), glue(v.1.8.0), parallelly(v.1.45.0), codetools(v.0.2-20), distributional(v.0.5.0), stringi(v.1.8.7), gtable(v.0.3.6), later(v.1.4.2), pillar(v.1.10.2), htmltools(v.0.5.8.1), brio(v.1.1.5), R6(v.2.6.1), textshaping(v.1.0.1), rprojroot(v.2.0.4), evaluate(v.1.0.4), shiny(v.1.11.0), lattice(v.0.22-7), R.methodsS3(v.1.8.2), bslib(v.0.9.0), httpuv(v.1.6.16), Rcpp(v.1.0.14), nlme(v.3.1-168), mgcv(v.1.9-3), xfun(v.0.52) and pkgconfig(v.2.0.3)